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SUMMARY 


This paper describes a simple numerical model for hurricane track prediction which uses a 
multigrid method to adapt the model resolution as the vortex moves. The model is based on 
the modified barotropic vorticity equation, discretized in space by conservative finite differences 
and in time by a Runge-Kutta scheme. A multigrid method is used to solve an elliptic problem 
for the streamfunction at each time step. Nonuniform resolution is obtained by superimposing 
uniform grids of different spatial extent; these grids move with the vortex as it moves. Preliminary 
numerical results indicate that the local mesh refinement allows accurate prediction of the 
hurricane track with substantially less computer time than required on a single uniform grid. 


INTRODUCTION 


Accurately predicting the track of a moving hurricane is a problem of great practical 
importance. One approach is to treat the problem as one in computational fluid dynamics, taking 
observed meteorological data as initial values for a numerical model. Many factors influence 
the accuracy of this approach, including the initial data (or lack thereof), the dynamical and 
physical processes included in the model, and the numerical scheme employed. While the relative 
importance of these three factors is a subject of considerable debate, in this paper we focus on the 
third. 

Our premise is that predicting the track of a moving hurricane accurately requires resolving 
the flow field adequately on both the large scale surrounding the vortex and the small scale within 
the vortex itself. Since the spatial scales involved may differ by more than an order of magnitude, 
models using uniform resolution are inherently less efficient than what should be possible. Here, 
we use a simple dynamical model which has been used successfully by many authors (ref. 1, 2, 3), 
namely, the modified barotropic vorticity equation. However, rather than use a single uniform 
grid as in those studies, we investigate the use of adaptive multigrid techniques, with the goal 
of obtaining high accuracy at low computational cost. In the following sections we detail the 
formulation of the model, describe the mesh refinement scheme, and present some preliminary 
numerical results. 
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MODEL FORMULATION 


Governing Equations 


We formulate the model on a section of the sphere using a Mercator projection (true at latitude 
$ — <p c ). The model consists of the modified barotropic vorticity equation 

^ + m 2 J(ip, C) + = um 2 V 2 (, (!) 

where the relative vorticity £ and streamfunction ip are related by 

(m 2 V 2 — 7 2 ) ip = C- (2) 

Here V 2 = d 2 /dx' + d~/dy 2 , J{ip, C) is the Jacobian of (ip, <) with respect to (x, y), /3 = 2Q cos (p/a 
(with a and O the radius and rotation rate of the earth), and m = cos <p c / cos (p is the map 
factor. There are two quasi-physical parameters: the diffusion coefficient v , and the parameter 7 
(inverse of the effective Rossby radius) which helps prevent retrogression of ultralong Rossby waves 
(ref. 4). We also consider versions on the /-plane (m = 1 and /? = 0) and /3-plane (m = 1 
and /? = 2D cos <p, /a). The model domain is a rectangle in x and y centered at (x, y) = (0, 0), 
where (A, (p) — (A, , <p ( ). At the boundaries we specify the streamfunction ip (and thus the normal 
component of the velocity); where there is inflow, we also specify the vorticity C- 


Space Discretization 


On a single uniform rectangular grid Q h consisting of gridpoints ( x,,y ,) with mesh spacing h in 
x and y, we discretize (1) and (2) in space by finite differences as 


+ ™]\j(ip, C) + Pjrrijdf'ip,., = i/m;V 2 ,0., 


(3) 


and 

(mjVjj - 7 ") fpKj = C,n ( 4 ) 

respectively. Here i.j(ip , C) the discrete Jacobian of Arakawa (ref. 5), and dj !l 'Uh t and 

V 2 ip tJ are the 0(h 2 ) centered difference approximations to dip/dx and the Laplacian operator, 
respectively. We apply (3) and (4) at the interior points. At boundary points where there is inflow, 
£ is specified; otherwise, we predict £ on the boundaries by applying an equation of the form (3) 
but using appropriate one-sided differences. It should be noted that using the Arakawa Jacobian 
is crucial here: the fact that it conserves discrete analogues of vorticity, enstrophy (mean square 
vorticity), and kinetic energy implies that the model is not subject to nonlinear computational 
instability. 
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To write the space-discretized equations in a more compact form, we collect the values ip,., and 
C , into grid functions ip h and , respectively, defined on the grid Q 1 ' . We can then write (3) and 
(4) as 

^ = F”(rP“, C") (5) 

and 

G h {ip h ) = C\ ( 6 ) 

where the operators F 1 ' and G 1 ' express the space discretization described above. 


Time Discretization 


To discretize (5) and (6) in time we use the classical fourth-order Runge-Kutta (RK4) scheme. 
To describe it, we specify a time step At > 0 and introduce time levels t* = kAt for k = 0, . . . - 
Suppressing the superscript h for simplicity, we now use the superscript k to denote values at time 
level k , e.g., ip 1 ' « ip h {h). With this notation, the RK4 scheme can be written as 


where 


C*' + ^ - C A = f a 


<1 

—•I'M 


£*'+5 — 

c* 

7 At 




At 


£< + l _ 

C A 


At 


:= C*), 

:= F(^ A+ 2,C A+ 5), 

:= F(/ + 5,C i+ ^). 


G(ip k+ *) = C A+ ^ 

G$* +1 ) = C* +1 . 

G{ip l+} ) = C k+ \ 


pk + 1 


1 ( F k + 2 F a+ 5 + 2 F k + ? + F k + l 

6 V 


(7) 


( 8 ) 


Thus, to execute a single time step t/, — > t* +1 , we perform the four stages indicated in (7); each of 
these stages consists of computing F based on known values of ip and £, predicting a new vorticity 
and solving the diagnostic equation for the corresponding streamfunction ip. 

Although it requires four times as much work (per time step) as the second-order Adams- 
Bashforth scheme commonly used in such models, this RK4 scheme has several advantages. First, 
it allows time steps at least four times as large, so in fact it is more efficient. Second, it is more 
accurate, so time discretization errors are less likely to distort the conservation properties of the 
Arakawa Jacobian. Finally, since it is a one-step scheme, it has no computational modes and needs 
no other method for the initial time step. 
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Multigrid Solution 


To solve the diagnostic equation at each stage for the streamfunction ip, we use a multigrid 
method. For the relaxation scheme we use a point Gauss-Seidel method formulated as follows. The 
discrete (interior) equation (4) can be written as 


(I4hi = p - Si.,) = -H = F,,, 


( 9 ) 


where 

S>.j := ipi - ij + Tpi+l.j + V'fJ-l + ' i Pi.j + i 

is the sum of the neighboring values of ip and 


a, := 4 + 


7 : V 


m~: 


( 10 ) 


( 11 ) 


is the diagonal term of the discrete Helmholtz operator. Given an approximate solution ip of (9), 
we relax at a point (i,j) by changing the value there to satisfy the corresponding equation (9); this 
results in the new value 


ifoj = 


VFij + S.j 


( 12 ) 


where S,j is defined using the current surrounding values in (10). The corresponding residual (if 
needed) is given by 

r,.j := F,.j - p {(Tjipi.j - S ^ - ip,.i) - (13) 

We use this relaxation (with red-black ordering) as a smoother in a multigrid method, using half- 
injection for the fine-to-coarse transfer of residuals and bilinear interpolation for the coarse-to-fine 
transfer of corrections. For the control algorithm we use repeated V(l,l)-cycles. 


LOCAL MESH REFINEMENT 


Given the premise that the flow near the center of the vortex requires much higher resolution 
than the flow surrounding the vortex, we now consider how to provide such variable resolution. 

Our basic method is essentially that of (ref. 6) , constructing nonuniform resolution by 
superimposing uniform grids of varying spatial extent. Since all calculations are carried out on the 
uniform grids, programming remains relatively easy. 

To illustrate the method, let us consider first the case of two grids: a coarse grid covering 
the whole domain fi, and a fine grid which covers only a portion of the domain (i.e., enclosing 
the vortex). We assume that the boundaries of the fine grid coincide with coarse grid lines. The 
model variables £ and ip are carried on both the coarse and fine grids (denoted by ( 2l> , ip 2h and ^ , 
ip h , respectively). Noting that the coarse grid allows time steps twice as large as those on the fine 
grid, we use the following basic procedure for stepping the model from time t * to t* + ] : 
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1 . Execute one time step of length At on the coarse grid to produce ; 1 , ip- 1 ' 1 -” 1 ; 

2. Execute two time steps of length At/2 on the fine grid to produce £ 

using boundary values for ip interpolated from the coarse grid (in space and time) ; 

3. Copy the fine-grid solution to the coarse grid at points common to both. 

Several points deserve mention here. First, in solving the implicit problem for ip on either grid, 
we use the multigrid method outlined above. This introduces additional coarse grids, e.g., a grid 
with mesh spacing 2 h covering only the region of the local fine grid S7 / '. In fact, the “underlying” 
part of the coarse grid f could be used for this; however, the resulting complications of 
preserving interface values (for fine-grid boundary values) and restricting relaxation to only part 
of Q J ‘ seem too high a price to pay for the relatively small savings in storage which would be 
achieved. Second, after completing the above three steps, the resulting solution on the composite 
grid Q- = fl h U ft 21 ' could be further refined by applying a composite-grid discretization of the 
governing equations; this FAC (Fast Adaptive Composite grid) method and several variants are 
described in (ref. 7), and will be explored in future work. Finally, the above approach generalizes 
immediately to more than two grids. 

For the initial work reported here, we have made the following simplifying assumptions. First, 
we require the grids to be rectangular and strictly nested (i.e., any fine grid is contained wholly 
within the interior of the next coarser grid), with one grid per level (i.e., the refinement occurs in 
one region only, surrounding the vortex). Second, we use a constant mesh ratio of two (i.e., the 
mesh spacing h on any grid is twice that of the next finer grid, if any). Finally, we will specify 
the number of grids and their sizes in advance but allow them to move following the vortex as the 
solution is computed. 

Since the problem to be solved has an easily identifiable region of interest surrounding the 
vortex, we take the following simple approach to moving the grids. First, we locate the vortex 
center on the finest grid. Then for each grid in turn, from the next-to-coarsest to the finest, we 
decide whether or not to move the grid. This decision is based on the distance of the vortex 
center from the center of the grid: if it is more than a specified fraction a of the distance L to the 
boundary, we move the grid. The move is calculated so as to “overshoot” a bit, i.e., aiming to put 
the vortex center beyond the (new) grid center by a specified fraction 6 of the distance to the grid 
boundary. Note that care must be taken at this stage to ensure the strict nesting of grids assumed 
above. Finally, the grid is moved by shifting the values which remain on the grid and filling in 
the rest by interpolation from the next coarser grid. For the results presented here, we check for 
possible grid moves after each time step on the coarsest grid, and use the parameters a = 0.4 and 
6 = 0 . 2 . 

To locate the vortex center (needed both for moving the grids as described above and for 
determining the vortex track), we first locate the point of maximum vorticity on the finest grid. We 
then interpolate the vorticity at that point and its nearest neighbors in x and y (five points total) 
by a quadratic function, and define the vortex center to be the location of the maximum of that 
quadratic. 
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RESULTS 


The initial conditions for the test problem consist of an axisymmetric vortex superimposed on 
an environmental flow, as considered in (ref. 1). The environmental flow is given by 

*■>-(¥)-(")• »■> 


which corresponds to the zonal current 


dip 

u{y ) = — j- = h (J sin 
dy 



(15) 


The tangential wind in the initial vortex is given by 


ne = 2K, (£) 


exp[-a(r/r„ l ) , ’] 

1 + ( r / r „,) 2 ’ 


(16) 


where r = [(x — Xu ) 2 + (y — yo ) 2 ] 1 ^ 2 is the radial distance from the vortex center (x n , yo)- Note that 
V has the approximate maximum value V m near r = r,„ (exact when a = 0); the exponential factor 
is included to make V vanish quickly for large r. The vorticity corresponding to (16) is 


C « = 


d(rV) 

rdr 


7 TOT -<~J 


(17) 


We will use the following parameter values: uo = 10 ms 1 and L — 4000 km for the environmental 
flow, and V m = 30ms ~\r m = 80 km, a = 10"°, and b = 6 for the initial vortex. The 
computational domain is a square of side length 4000 km on a /5-plane, using /5 for the latitude 
20° N; the vortex is initially centered at x 0 = 750 km and y„ = -750 km. The model was run from 
t = 0 to t = 72 hr; for simplicity we have set v — 0 and 7 = 0 here. 

To establish a standard for comparison, we ran the model with high resolution (384 x 384 grid 
with spacing h = 10.42 km and time step 10 s). We then ran the model with a variety of grid 
configurations (using up to four grids) and compared the vortex track to that of the reference run. 
Table I summarizes these results, with the runs listed in order of increasing execution time (on a 
SUN SPARCstation2) . All of the cases in this table use only square grids, with N, = N y = N. The 
forecast error is defined as the distance between the predicted vortex location at a given time and 
that in the reference run. These results show that the local refinement process has the potential to 
substantially reduce the execution time required to achieve a given accuracy. For example, a sing e 
grid with h = 31.25 km (run 6) achieves errors on the order of 10-20 km; with local refinement 
(run 2) comparable accuracy is obtained with only about 36% as much computer time. Similarly, a 
single grid with h = 20.83 km (run 8) achieves errors on the order of 1-5 km; with local refinement 
(run 7) comparable accuracy is obtained with only about 42% as much computer time. In fact, ^ 
run 7 with local refinement achieved about the same accuracy as did the single-grid run with h = 
15.625 km (run 9) but with only about 18% of the computer time. In addition, the solution fields 
produced with local refinement (run 7) are smooth, as shown in Figures 1-5, with no indication of 
any problem due to the change of resolution at the grid interfaces. 
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CONCLUDING REMARKS 


The preliminary results reported here show that adaptive multigrid techniques can substantially 
reduce the computer time required to make accurate hurricane track forecasts. In addition to 
ongoing testing of the existing model, we plan to investigate the following possible improvements. 
First, we plan to include the FAC method as discussed above. This should have the advantage of 
more precise conservation of vorticity, enstrophy, and kinetic energy at the grid interfaces. Second, 
we plan to construct a fully adaptive version of the model by using the Full Approximation Scheme 
(FAS) to produce estimates of the local truncation error to be used in an automatic grid refinement 
scheme (as proposed in ref. 8) . Finally, we plan to test the model using real data, and compare its 
performance to that of models currently in operational use. 
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Table I. Results of Model Runs 



Grid 

size{s) 

At 

Forecast error (ki 

m) at: 

Execution 

Run 

N 

h (km) 

(min) 

24 hr 

48 hr 

72 hr 

time (sec) 

1 

64 

62.50 

60 

110 

143 

47 

170 

2 

64 

62.50 

60 

11 

8 

17 

504 


64 

31.25 

30 





3 

96 

41.67 

30 

53 

12 

25 

799 

4 

32 

125.0 

120 

14 

24 

39 

916 


32 

62.50 

60 






48 

31.25 

30 






64 

15.62 

15 





5 

64 

62.50 

60 

1 

6 

10 

1,174 


64 

31.25 

30 






64 

15.62 

15 





6 

128 

31.25 

30 

11 

8 

19 

1,409 

7 

64 

62.50 

60 

1 

5 

5 

2,047 


64 

31.25 

30 






96 

15.62 

15 


— 



8 

192 

20.83 

20 

1 

3 

5 

4,860 

9 

256 

15.62 

15 

2 

3 

4 

11,405 

10 

384 

10.42 

10 

- 

- 

- 

41,716 
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